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ABSTRACT 



A comprehensive computer code has been developed through the use of the 
governing equations of motion and the Boussinesq hypothesis and the rise of a vortex 
pair in a linearly density-stratified medium has been calculated through the use of the 
Biot-Savart law and a finite difference scheme. The results have shown that the code is 
capable of predicting the behavior of the vortices and the propagation of the internal 
waves resulting from the vortex motion. The reasons for the finite rise of the vortices 
in a stratified medium is clearly explained in terms of the counter-rotating vortices 
induced by the stratification. 
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I. INTRODUCTION 



Vortices and vortex wakes have become a major theme of aerodynamics research 
since the advent of the large aircraft and the understanding of their evolution required 
an examination of many fundamental problems in .fluid mechanics. Much of the 
progress made during the past two decades was discussed at the Symposium on 
Aircraft Wake Turbulence and Its Detection [Ref 1] and at the Aircraft Wake Vortices 
Conference [Ref 2]. Comprehensive reviews of the entire subject have been given by 
Donaldson and Bilanin [Ref 3], Widnall [Ref 4] and Hallock and Eberle [Ref 5]. 

These studies, as well as numerous others carried out since 1977, have uncovered 
a number of complex problems which must be resolved in order to achieve a better 
understanding of the important features of -trailing vortices in homogeneous and 
stratified media. The principal ones are as follows [Ref 6]: 

(a) Roll-up process; The velocity and turbulence distribution at anv station 
behind the wing depend on the wing section, wins-tip shape^. Revnolds 
number, wins incidence, and the distance of the stationTrom the wing [Ref 7]. 
The distribufions of the initial velocitv and turbulence, which influence the 
roll-up and the decav process, cannot be changed independently. For 
example, a chanse in tip shape changes the core size, as well .as the velocity 
and turbulence distributions. Hish levels of turbulence result m an increased 
diffusion of vorticity, which in turn increases the core size. 

(b) Probe sensitivity of the vortices: Flow visualization studies suggest that 
trailing vortices 'are extremelv sensitive to disturbances created .by "even veiw 
small probes or bubbles. This forces one to use non-intrusive means . of 
measurements such as an LDV. Even then, 'vortex wandering' [Ref 8], which 
makes the vortices appear larger than normal in time-averaged velocity 
measurements (for vortices senerated by a wing in a wind tunnel), or the 
unsteadv nature of the flow ffor vortices' senerated bv a wing in a tow basin) 
makes the mean velocity prohles in the vortices difficult to determine. 

(c) Large-scale instabilities: The vortices are seldom observed to decay aw-ay 

• owihs to viscous and turbulent dissipation, but are almost always. 'destroyed 

bv either mutual induction instability (Crow instability [Ref 9]) and/or vortex 
breakdown. The Crow instability grows exponentiallv, and results either in a 
linking of the vortex pair into a series of crude vortex rings or in a hishly 
disorsanized interminslins of the vortices. Vortex breakdown, wlios'e 
mathematical details have" not yet been adeouatelv treated, rearranges the 
vortex structure and increases the' core size, turbulence, and energy dissipation. 
Thus, it is verv difficult to rneasure accurately th.e trajectories"bf the three- 
dimensional vo'rtices from their creation to their ultimate demise. 

(d) Reynolds number; Even the highest Revnolds numbers, based on wimi chord, 
rea'ched in wind tunnels or tovinsz basins, are in order of maanatude lower 
than what is possible for an aircraft. Thus, the scale effects are not easy to 
assess. 

(e) Ambient conditions such as turbulence and stratification play major roles in 
the evolution of vortices. The quantification of these 'effects requires 
numerical analysis and extremely careful experiments. 
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(f) Ground or free-surface efTects: The vortex, pair mav move towards a risid 
boundary' at which, the no-slip condition must be satisfied or towards a free 
surface at which the zero-shear condition must be satisfied. In either case, the 
vortices come under the influence of their images and move accordingly. 

The phenomenon is further complicated by several additional facts. When the 
vortices are propelled towards a rigid surface, vorticity of opposite sign is generated on 
the no-slip boundary* and swept towards the vortex pair. The total vorticity diminishes 
very quickly as vorticity from the two regions diffuses, the wall region serving as a 
strong sink for the vorticity associated with the -original vortex [Ref 10]. The 
developement of a boundary layer along the rigid wall may give rise to flow separation 
for sufficiently high Reynolds numbers. With or without such a separation, however, 
the center of the vortex pair eventually moves away, or' 'rebounds', from the wall 
[Refs. 10-12]. 

For the case of a zero-stress boundary, the free surface still acts as a vorticity 
sink, but this is relatively weak due to the absence of intense oppositely-signed 
vorticity. Thus, in the absence of other impeding phenomena, one expects a mild 
interaction betw'een the vortices and the free surface and a small rebound of the vortex 
center from the free surface. However, the ability of the free surface to deform under' 
the influence of strain fields leads to a strong interaction between the vortices and the 
free surface. 

It is evident from the foregoing that the motion and the life-span of trailing, 
vortices are goverened by a number of nonlinearly-dependent complex phenomena.' A 
number of experimental and analytical studies have been carried out at the Naval 
Postgraduate School by Sarpkaya and his students [Refs. 6,13-17] in order to 
investigate the efTects of these parameters on the rise and demise of the trailing vortices 
in homogeneous and density stratified media. These studies have clearly identified the 
various demise mechanisms in both media and established basic relationships between 
.the rise of vortices and the governing parameters in a finite as well as effectively infinite 
medium [Ref IS], free from ambient turbulence. The present investigation is a 
continuation of the foregoing studies and is confined to the computer modeling, via a 
finite differencing technique, of the rise and demise of trailing vortices in stratified 
media. 
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II. EXPERIMENTAL EQUIPMENT AND PROCEDURES 



A. EQUIPMENT 

The equipment used to generate the trailing vortices has been extensively used at 
this facility over the past three years [Refs. 6,14-19]. Only the salient features, most 
recent modifications and the adaptation for this work are briefly described in the 
following. 

The system consists essentially of a towing basin. The auxiliary components 
consist of the brine fill system, turbulence management system, top and bottom 
carriages, velocity measuring system, hghting system, and the model [Refs. 14-19]. The 
brine fill system is controlled by a computer to provide a consistently repeatable 
density stratified media [Ref. 19]. 

Drains are provided at the bottom of each end of the basin. Two parallel rails 
are mounted along the bottom of the tank. A carriage rides smoothly on these rails 
and provides the test body with a constant velocity through the use of an endless cable 
and a DC motor. The velocity of the model is measured and continuously monitored 
through' the use of a magnetic linear displacement transducer. It yields a signal 
proportional to the velocity of the model with an error of less than one percent. 

The two rails, the carriage and the filling pipes are located on or near the bottom 
of the basin and under a turbulence management system. This system consists of one 
inch thick polyurethane foam sandwiched between two perforated aluminum plates. 

B. MODEL 

A rectangular foil (NACA 0012) was used in the present study. The fojl had the 
following dimensions: B = 4.50 in., b^ = 3.53 in., and c = 2.32 in. The interior of 
this model vvas hollowed and used as dye reservoir to seed the vortex cores. The model 
was mounted on its base by means of a thin streamlined aluminum bar with a cross 
section of a NACA 0006 foil and set at the desired angle of attack. As noted earlier, 
the model was pulled by means of a DC motor, pulley, and cable system at the desired 
speed. Additional details of the model construction and mounting are discussed by 
Johnson [Ref. 15]. 



C TEST PROCEDURES 



A model is placed in the basin and the basin is filled gradually with a brine 
solution [Ref. 19] to the desired le\-el. In the cases where dye was used to seed the 
vortex pair, the model was filled with dye prior to filling the basin. .After removal of 
trapped air and after sufficient time to allow for the escape of dissolved air in the water 
and the elimination of any internal currents in the basin, the model was set in motion. 
Most of the experiments were repeated at least three times. 

The motion of the trailing vortices were recorded on highspeed Polaroid film at 
the test section (one of the plexiglass panels near the middle of the basin). Each 
picture included tw'o clocks accurate to 0.1 seconds, the vertical and horizontal scales 
on the plexiglass window and, of course, the side view of the trailing vortices as they 
rose from the model after formation. The time interval between successive pictures is 
determined from the two clocks. The first picture always captured the model the 
instant it passed through the test section. The subsequent pictures (taken at about 
0.75 second intervals) captured the rise and demise of the vortices. The vertical rise of 
the vortices is determined from the verticle scale. Attention has been paid to the fact 
that the vortices are farther away from the scale on the window and that the scale 
placed vertically in the middle of the test section does not exactly correspond to the 
scale marked on the window due to refraction and parallax. The necessary correction 
was made by photographing a scale placed in "water in the middle of the test section 
together with a scale marked on the window*. This resulted in a simple conversion 
table which enabled the determination of the actual position of the vortex from the 
scale reading on the photograph. 

The results are normalized and plotted in various forms and compared with those 
obtained in the previous runs. Each experiment w*as repeated three times for most of 
the data presented herein. All trailing vortices w'ere recorded on film until the time 
they have completely dissipated either due to aging (diffusion of vorticity due to 
viscousity, turbulence, entrainment, and detrainment), or due to Crow instability 
(sinusoidal instability and linking of vortices), and, or due to vortex breakdown (core 
bursting). It was thus" possible to determine the life span of vortex cores from their 
formation to their final demise. 
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D. COMPUTER MODEL 

The computer model (Appendix A) is run on the IBM 3300 main-frame computer 
at the Naval Postgraduate School. The computer code was written to run on either the 
MVS Batch Processing system or the VM,CMS system. The stratification parameter 
specified for the experimental runs is used to determine other parameters [Ref 1'] used 
by the program. The code is compiled under VS Fortran and is loaded and executed 
using DISSPLA version 10.0. The code was \^Titten such that the user has the capacity 
to stipulate the number of time step iterations that are executed. .Also, the time steps 
at which a constant density contour, complex velocity profile, density perturbation 
contour and the vorticity contour are plotted,- are designated by the user prior to 
compilation. These plots are analized to determine the normalized rise height of the 
initial vortex pair until its demise. These are then compared with those obtained 
experimentally. 

The computer model is an intense computational algorithm that requires 
approxiamately 12.5 minutes of Central Processing Unit (CPU) time per iteration. Due 
to this high use of CPU time, the code has been designed to he executed until stoped at 
a user defined step iteration and then save all necessary data required to restart the 
model. The code can then be restarted at that step iteration at a later time. The 
program VORRS (VORtex ReStart) has been specifically w'ritten to accomplish this 
task. 
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III. ANALYSIS 



A. NUMERICAL ANALYSIS OF THE VORTEX .MOTIONS 

The seneration of the internal waves and the rise and demise of the vortices in a 
stratified medium may well be analyzed through the use of the equations of motion for 
an incompressible fluid for both laminar and turbulent motions provided that a suitable 
turbulence closure model is adopted and the usual Boussinesq approximation 
(gravitational acceleration is much larger than the fluid accelerations) is made. For the 
type of motions considered herein the Boussinesq approximation is quite valid and has 
been used in the investigation of all types of internal waves in stratified fluids. 

For a two-dimensional flow, with y vertical and x horizontal, the equations of 
motion are 



5u du 5u 1 5p • T 

— + u — + V — = + vV^u 

5t dx dv p dx 



and 



dv d\ dv 1 5p , 

h u — + V — = - g + vV'^v 

5t dx dy p dy 



( 3 . 1 ) 



( 3 . 2 ) 



and the equation of continuity 



5u dv 

— + = 0 . 

dx dy 



( 3 . 3 ) 



Defining vorticity as usual by 



3u dv 

dy dx 



(3.4) 



and eliminating pressure between Eq. 3.1 and Eq. 3.2 , one has 
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(3.5) 



_L H Z H Z 

di dx dy 



vV^ ^ 



_g__£p' 

Po 



in which p' is the fluctuating part of the density p given by 



P = Po P(y) ^ P'(x,y,t) 



(3.6) 



where p^ is the reference density, and p(y) the initial density at elevation y at t = 0. 
The diffusion of density is given by 



d d d 

( — + u — + V — 
5t dx dy 



- vV^)p = 0 



(3.7) 



Combining Eqs. 3.3. 3.6 and 3.7, and simplifying, one has 



dp' dup' dyp' 

— ^ 4- - 

di dx dy 



V 



— + v7-p' + V' 
dy 



d^p 
__ _ 

t7y- 



(3.8) 



It is convenient to cast the foregoing equations into nondimensional forms, 
scaling each variable by a quantity characteristic of its expected magnitude. Two 
possible time scales exist. There is the dynamic time scale, which is the time a 
characteristic length would be traversed by a fluid particle traveling at the characteristic 
velocity, and there is the buoyant time scale based on the natural buoyancy frequency 
of the stratified flow, i.e., the Brunt-Vaisala frequency N, as defined by Daly [Ref 17]. 
Each of the scales gives a slightly different form of the normalized qoverning equations. 

Dynamic Scale: 

Introducing and as the characteristic velocity and length, one has 



^L, U t U L U 

Cm = tm = — , Re = F = — 

L ’ V ’ (eL)' 



Normalizing Eqs. 3.5 and 3.8 



m 



5Cm ^^mCm 5V C 1 . 1 C* P 

- -m ^ ( m -m ^ m -m \ _ ^2^ j _ 

5t ^Xm Re F - dx 

ni m ^ rn v m 



(3.9) 



(3.10) 
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and 



dp ' p_' dvo' ^ 1 , d^p^ 

Jim ^ ^ m£m j ^ ^ -—S'-) 






m 



dx 



m 



0 \' 



m 



Re 



m 



OX 



m 



(3.11) 



in which 



^ U ’ 

F. .- = — ^. Re = -^. n = ^ 



sL. 



L>. 

V 



N 

( 

N 



a 

o 



N/ = — , N 2 = 

C T ’ O 



g p(y) 
p. ■ 



g ^Pm _ 



L dy 



m 



Equations 3.10 ans 3.11 are valid when > > 1, i.e., when the bouyancy has 
ver>' little influence on the nonlinear dynamics of the motion. In this case, the density 
perturbation acts as a passive scalar advected by the velocity field. 

Bouyant scale: 

Introducing the following dimensionless parameters; 



V = 
' m 



Pm - ’ TT * TT ' 



m 



= 



u 



m 



V 

u" 



m 



X 

r 



'^'c. ‘Sm 



C ■■ 



T 2 _ 



L. 



(3.12) 



and normalizing Eqs. 3.5 and 3.8, one has 



dC du C dv C 

~m p ^ m •‘m m ^ 



di 



m 



d\ 



in 



’ dx 



= i + -^1 



m 



Re 



cx 



m 



(3.13) 



and 



dp^' ^u_p_ f5v p 

j_ 17 / m'^m ^ _ p ^2.^ 



<2t 



+ Fv(- 



m 



dx 



m 



dx 



m 



m 



(3.14) 



in which F^,, Re and n" are as defined in Eq. 3. 1 1. 
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Equations 3.13 and 3.14 are valid when < < 1 and buoyancy dominates the 

flow. When approaches zero,, the equations that result from Eqs-3.13 and 3.14 

describe the propogation of linear internal waves. The F^. < < 1 regime is of interest 
in the present investigation because for submerged bodies of naval interest F^. (the 
Froude number) is about 0.01. The analysis will consider the full nonlinear equations 
given by Eqs. 3.13 and 3.14 rather than their linearized form (i.e.. the case of F^. = 0). 
However, the flow will be assumed to be inviscid, i.e., the viscous difiusion will be 
ignored to a first order approximation. Subsequent analysis will incorporate the eflects 
of viscous and turbulent diffusion into the numerical calculations. The velocities u and 
V are given by the Biot-Savart law as 



u(x,y) 



J 



(y^ - yK(x\y')dx'dy^ 
2jtr^ 



(3.15) 



and 



u(x,y) = J 



.(x - x')C('x',y')dx'dy' 
27tr^ 



(3.16) 



in which (x,y) is the point at which the velocity is calculated and (x',y') is an arbitrary’ 
point at which the vorticity is ^(x',y'). 

The radial distance r is given by: 



r^ = (X - X')- + (y - y')- . 



(3.17) 



Normalizing Eqs. 3.15 - 3.17 through the use of the characteristic dimensions 
given in Eq. 3.12 






27ir * 
m 



(3. IS) 



F u = r 

V m J 



27ir - 

m 



(3.19) 



and 
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(3.20) 



It is appropriate to take V^, the initial mutual induction velocity of the vortex 
pair, as the characteristic velocity and the b^, the initial separation of the vorticies. as 
the characteristic length. Thus, one has: 



L 



c 






V 

Fv = ^ 



X b 

C 0 



V 

o 




In dimensionless forms, one has: 



(3.21) 



N b 



O _ 



= 27tF. 



n = F ^ ° 



c o 



V 



t^ = tN = — ^ 
m b F 

o V 



(3.22) 



in which N^b^/V^ is called the stratification parameter SP. The value of SP normally 
varies from zero to about 100. The larger the SP, the stronger the stratification. For 
example, for a model running at a speed of 5 ft 'sec at an angle of attack (a) of 10°, a 
SP‘ smaller than 1.0 represents weak stratification, a SP between 1 and 5 represents 
medium stratification, and a SP larger than 5 corresponds to strong stratification. 

B. COMPUTER AIDED MODELING OF VORTEX MOTION 

The computer algorithm for the numerical solution of Eqs. 3.13 and 3.14 is 
provided in Apperidix A. The solution is done assuming the case of inviscid flow and 
the initial vorticity distribution was assumed to be Gaussian, i.e.: 




F ' r ^ • 
-^exp(- :^) 



(3.23) 



A computational domain of 5b^ by lOb^ is used in the algorithm in order to 
mimick the physical features of the tow basin. The’ no-slip condition is satisfied on the 
walls corresponding to the chanel walls. The zero normal velocity condition is assumed 
at the free surface and along the vertical line of flow symmetiy. All computations are 
done in the fourth quadrant. The fields of velocity, vorticity. density, and the 
fluctuating components of the density are calculated at each time step and plotted at 
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regular intervals. The computer algorithm is quite robust and can be applied to any 
stratified and unstratified flow situation. As the velocity field is computed from 
analytic expressions (i.e., the' Biot-Savart equations) vice an algorithm involving an 
iterative scheme, stability of the program is enhanced at a slight expense in execution 
time. An upwind finite differencing technique [Ref 20] was used with an original grid 
dimension of 24 by 44 nodes with a spacing of0.2ob^. This grid size was subsequently 
increased to 34 by 64 and 44 by 84 in order to examine the dependency of the results 
on the grid size. 

The mesh coordinates (x,y) covering the computational domain (of dimension 
m-4 by n-4) are: 



X e 2,3,4,...,m-2 

(3.24) 

Ye {2,3,4,.. .,n-2). 

The necessity for the the overlap of two nodes on all sides of the computational 
domain is to facilitate satisfying the boundary condiditons. Eqs. 3.13 and 3.14 are 
solved over the original mesh which results in vorticity and density information. 

The two complex velocity grids are derived by shifting the original mesh h/2 in 
the X direction for the u velocity and h/2 in the y direction for the v velocity. The 
values of these new nodes are determined -using the vorticity derived' for the original 
computational grid using the Biot-Savart law, Eqs. 3.15 and 3.16. The velocities at the 
shifted grid points are then averaged across the computational grid points for further 
use in the finite differencing scheme. 
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IV. DISCUSSION OF RESULTS 



The computer code described in the previous section has been used to carry out 
calculations for four specific cases. The first of these dealt with the case of the rise of 
vortices in a homogeneous unstratified medium. The grid size as well as the core 
radius were varied within limits to e.xplore the differences in the motion of the vortices. 
Figures B.l through B.13 show the velocity and vorticity contours at various 
normalized limes. Clearly, the initial vonex spacing increases with time and the 
vortices remain symmetrical due to the conditions imposed on the numerical analysis. 
In other words, these vortices cannot be subjected to sinusoidal instability or core 
bursting. Thus, the numerical results should agree more closely with the experimental 
cases where the vortices did not undergo such instabilities. Figure B.14 shows the 
e.xperimental and numerical results. The experimental data are separated as 'HIGH' 
and 'LOW' because of the fact that the crests and troughs of the sinusoidal instability 
were separately evaluated. The numerical results for three normalized core sizes are 
also shown in Figure B.14. It is clear that the core size has a pronounced effect on the 
numerical predictions and that the results obtained with r^ = 0.08 agree more closely 
with the experimental data. As noted earlier, the numerical res-ults should be compared 
with the average instantaneous height rather than with the 'HIGH' and 'LOW' of the 
vortices. There is no conclusive analysis or data in the literature regarding the vortex 
core radius. It has been shown previously [Ref 16] that the core radius depends 
strongly on the shape of the trailing edges of the lifting surface. The better rounded 
the edges of the control surface the tighter the winding of the vortex sheets , and thus 
the smaller the vortex core. It appears that a vortex core of 0.07 or 0.08 should be 
used in future calculations. 

Figures B.15 through B.31 show lines of constant density, velocity, density 
perturbations, and vorticity contours for the statification parameter of 1.0. F^. = 
0.02976, r^ = 0.09, Y/b„ = -8.0 and X.b = 0.5. The srid used was 44 bv 84. 
These figures show that at small times the vortices rise vertically and the circulation in 
the flow field is primarily due to the initial vortices. As time increases (see e.g., Figure 
B.21b), the regions of circulation with countersigned vorticity develope in the upper 
right and lower left regions of the vortex. It is easy to show that this counter vorticity 
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not only reduces the rise velocity of the original vortex pair but pushed the pair against 
each other. Consequently, the countersigned vorticity begins to dominate the flow (see 
e.g., Figure B.2Sb) and the vortex pair migration stops. With further increases in time 
the vortex pair begins to migrate downward. 

The density contours reveal the same phenomenon in a difTerent context. .A,s the 
vortices rise, fluid of greater density is pushed upwards (see e.g., Figure B.25a) into 
regions of lesser density. Since such a migration cannot go on indefinitely, the vortices 
rise to a maximum height and then begin to sink downwards. The calculations do not 
take into account the sinusoidal instability and the vortex breakdown. Consequently, 
the vortices in the numerical calculations continue to exist unimpeded by the instability 
mechanisms. In reality, the vortices begin to break up as they near the end of their 
maximum migration and eventually disappear. 

The experimental and calculated values are compared in Figure B.32. The 
correspondence between the measured and calculated values is surprisingly good up to 
the time of maximum rise. This is partly because of the experimental fact that the rise 
of vortices in a highly stratified medium (N b /V = 1.00, F„ = 0.02976) is not 
strongly dominated by sinusoidal instability or vortex breakdown. The migration of 
vortices is inhibited primarily due to the reduction of vorticity of the initial vortices and 
the creation of countersigned vorticity. 

The numerical calculations were also applied to a medium with a relatively 
srnaller stratification for which the rise of vortices are closer .to that of the 
homogeneous case and thus the effect of the core radius is more pronounced. The 
results obtained with N' b^'V = 0.50, F„ = 0.02976 and with r = O.OS are shown in 
Figures B.33 through B.43. The observations made earlier regarding the velocity, 
vorticity, and density contours may be made with this particular stratification also. 
Obviously, the vortices rise to a greater height and the countersigned vorticity is 
relatively weaker and developes at a later stage. The calculations with this particular 
stratification have been repeated by varying only the initial core radius. Figure B.44 
shows the representative experimental data and the numerical results obtained with 
initial core- radii of 0.07, O.OS. and 0.09. Evidently and as anticipated on the basis of 
the previous calculations the initial core radius of 0.07 yields results which are in 
excellent agreement with those obtained experimentaly. Figure B.44 also shows that 
the rise and demise of vortices must be a strong function of their initial core radii. 
Experiments in both the laboratory and the field have shown that to be true and 
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pointed out the fact that the life of a vortex is dictated partly by the conditions 
prevailing at its creation and partly by the conditions surrounding the vortex at later 
stages of its motion. 

The results presented so far have dealt with vortex pairs released at relatively 
large depths of submergence and have shown that the computer code predicts results 
which are in good agreement with those obtained experimentally provided that a core 
radius of about r = 0.07 is used in the calculations. It has further been shown that 

O 

the smaller the stratification parameter the more dependent the results become on the 
initial core readius. Thus, the numerical analysis has not only provided a power of 
prediction for the motion of the vortices but also helped to delineate the proper value 
of the initial core radius. 

Having assertained the characteristics of the model and the features of the 
velocity, vorticity, and density variations it was deemed necessary' to explore problems 
of greater practicle and operational concern. It is because of this reason that the depth 
at which the vortices are generated was chosen to be Y/b^ = — 2.0 and the 

subsequent motion of the vortices was explored. Figures B.45 through B.52 show' the 
densitv, velocitv, densitv fluctuations and vorticitv contours at suitable time intervals. 
As the vortices approach the free surface and begin to move sidew-^ays under the 
influence of their images, the constant density contours begin to intersect with the free 
surface (see e.g.. Figure B.49a) and move sideways also. The intersection of these 
density contours is expected to give rise to intemal-wave-generated surface scars. The 
prediction of the characteristics of- these scars in terms of the vortices is of importance 
and requires that the free surface condition be changed. This change is to allow' for 
free surface deformation under the influence of the tangential as well as normal 
velocities induced at the free surface by the vortices. 

Additional calculations are underway to examine the effects of free surface 
deformation, viscous diffusion and ambient turbulence for various degrees- of initial 
stratification. The results presented herein and the initial comparisons are extremely 
encouraging and are expected to lead to a better understanding of this extremely 
complex and challenging phenomenon. 



V. CONCLUSIONS 



The investigation reported herein warranted the following conclusions: 

(1) The motion of trailing vortices generated by submerged bodies in a stratified 
or unstratified medium can be effectively modeled by the developed computer 
algorithm: 

(2) The computer model has aided in the delineation of the proper value for the 
initial vortex core radius; 

(3) The initial vortex core radius and vorticity strength are the primary factors 
governing the stability of the numerical results; 

(4) The migration of vortices is inhibited primarily due to the reduction of the 
vorticity of the initial vortices and the creation of countersigned vorticity; 

(5) The computer model could be made more realistic through the inclusion of a 
deformable free surface in the computations. 
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APPENDIX A 

VORTEX MOTION COMPUTER MODEL 



c 

C PR0GRA2-I VORTEX 

C 

c 

C 20 NOVEMBER 86 PROGRAMMER : BRIAN S. MILLER 

C 

c 

C PURPOSE : 

C 

C THIS PROGRAM NUMERICALLY SIMULATES THE RISE OF A TRAILING 

C VORTEX PAIR SHED OFF A LIFTING SURFACE SUBMERGED TN A 

C DENSITY STRATIFIED MEDIUM. 

C - 

c . . 

c TO RUN THIS PROGRAM: 

C 

C 1) GO TO "INPUT PARAMETERS" SECTION OF CODE AND CHANGE ANY PROBLEM 

C PARAMETERS ( SP , FV, DT, ZTMAX, RO , MP , NP , ISAVE) AS DESIRED. 

C 

C .2) SET MAXIMUM TIME IN DATA STATEMENT MAXSTP ; 

C EX: DATA MAXSTP /N/ 

e ( WHERE N= MAXIMUM NUMBER OF TIME STEPS) 

C 

C 3) SET INTERVALS FOR PLOTTING BY ENTRIES IN DATA STATEMENT NPLOT ; 

C EX: DATA NPLOT/ 1,2,8,12,20,40,50,75/ 

C ' ( THIS PLOTS OUTPUT AT THE 1ST , 2ND , 8TH , 12TH , 40TH , 50TH , AND 

C 75TH PROGRAM TIME STEP) 

C 

C 4) CHOOSE WHAT OUTPUT DEVICE TO USE, EITHER TEK 618 OR COMPRES 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



(METAFILE OUTPUT) BY GOING’ TO GRAPHIC DEVICE SECTION OF PROGRAI-I 
AND COMMENTING OUT WHICH OF TWO DEVICES YOU DONT WANT TO USE. 

5) COMPILE PROGRAI^ UNDER VS FORTRAN USING THIS COMMAND; 

FORTVS VORTEX (OPT (3) 

( NOTE: THIS OPTIMIZES CODE AT LEVEL 3 ) 

6) GO TO APPROPRIATE TERMINAL . AI'ID RUN UNDER DISSPLA. 

A. WHEN THE DISSPLA EXEC PROMPTS YOU ASK FOR 5 CYLINDERS OF 
TEMPORARY DISK STORAGE AND ALSO YOU MAY WISH TO ISSUE 

B. WHEN THE EXEC PROMPTS YOU YOU MAY WISH TO ISSUE: 

FILEDEF 06 DISK VORIOUT LISTING (PERM 
IN ORDER TO REROUTE THE PRINTED OUTPUT TO A DISK FILE 
NAMED: VORIOUT LISTING A1 . 

C. PROGRAM WILL NOW RUN SENDING OUTPUT TO YOUR DESIGNATED 
DEVICES. 



ADDITIONAL PROGRAM NOTES: 

1) THIS PROGRAM COMPUTES DENSITY AND VORTICITY IN A (M-2 X N-2) GRID 

BY SOLVING THE BOUYANTLY SCALED FORM OF THE NAVIER- STOKES 
EQUATIONS VIA AN "UPWIND DIFFERENCING" FINITE DIFFERENCE METHOD. 

2) THE VELOCITY FIELD IS COMPUTED USING A FORM OF THE BIOT-SAVART 

LAW AT EACH NODAL POINT INCLUDING THE CONTRIBUTION OF THE 
"IMAGE VORTICIES " ADDED TO INCLUDE THE FREE SURFACE EFFECT. 

3) THIS VERSION OF THE PROGRAM INCLUDES A DISSPLA ROUTINE 

TO PLOT THE DENSITY PERTUBATIONS , CONSTANT DENSITY, 

• VORTICITY, AND VELOCITY FIELD AT EACH SELECTED TIME STEP. 

THESE SELECTIONS FOR PLOTTING ARE INDICATED BY ENTRIES 
IN THE ARRAY NSTE? . MAXIMUM NUMBER OF TIME STEPS IS INDICATED 
BY AM ENTRY IN THE ARRAY MAXSTEP. 
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c 



c 

C 4) 

C 

C 

C 

C 

C 

c 

c 

C 5) 

C 

C 

C 6) 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 7) 

c 

c 

c 

c 

c 

c 

c 

c 

c 



THE COMPUTATIONAL GRID (M-2)X(N-2) IS USED FOR ALL COMPUTATIONS 
EXCEPT FOR VELOCITY CALLS TO THE 3IOT-SAVART SUBROUTINES. 

TWO ADDITIONAL GRIDS, ONE FOR U AMD ONE FOR V, ARE SHIFTED H/2 
UNITS IN THE X AND Y DIRECTION RESPECTIVELY. THE SUBROUTINES 
COMPUTE U AND V VELOCITIES ON THESE GRIDS USING VORTICITY FROM 
COMPUTATIONAL GRID. THE VELOCITIES ARE THEM AVERAGED ACROSS THE 
COMPUTATION.AL GRID POINTS TO BE USED FOR THE FINITE DIFFERENCING. 

THIS VERSION OF THIS PROGRAM WILL BE WRITTEN IN SINGLE PRECISION 
AND IN AN OPTIMIZED FORM OF THE ORIGINAL CODE. 

THIS PROGRAM HAS THE OPTION OF SAVING THE RHO AND ZETA ARRAYS 
AND OTHER PROGRAI-I CONSTANTS IN ORDER TO RESTART THIS SIMULATION 
AT THE TIME A PREVIOUS RUN HAD FINISHED. PROGRAM VORRS IS 
SPECIFICALLY WRITTEN TO INPUT THE DATAFILE WRITTEN BY THIS 
PROGRAM AND RESTART THE SIMULATION WITH THE SAME PARAMETERS. 

TO ACCESS THIS SAVE FEATURE SET VARIABLE ISAVE =1 . A FILEDEF 
OF FORM . . . FILEDEF 08 DISK <FN> <FT> <FM> . . . MUST BE ISSUED 
AT THE APPROPRIATE TIME IN EXECUTION. 

(NOTE: TO DISABLE THIS SAVE FEATURE SET ISAVE =0 ) 



TO CHANGE THE MESH DENSITY OF THE COMPUTATIONAL. DOMAIN THE 
USER SHOULD DO THE FOLLOWING: 

A. ADJUST M AND M 

B. ADJUST INITIAL VORTEX POSITION: MP , MP 

C. ADJUST THE DIMENSIONS OF THE W ARRAY IN SUBROUTINES 

CPLOT, ZPLOT, AND DRPLOT TO BE (MW, MW) 

D. IF THE GRID IS LARGER THAN 50 X 100 THEN ADJUST THE 

DIMENSION STATEMENTS IN THE MAIN PROGRAM AND IN 
SUBROUTINES ZPL, VPLOT AMD V2PLOT TO THE NEW GRID 
DIMENSION. 
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c 

c 

c 

Qr::k7:'k'k:k'k7^:k'k7:'k7^'k7^'k7^:k'k-k'k7^7:7^7z:k‘k'k7:'k'k'k'k'k:k7^:k:k:k'k:k:k'k:k'k'k'k'k^ 

c 

C MAJOR VARIABLES USED IN THIS PROGRAM; 

C 

C REAL VARIABLES; 



C T = 


NONDIMENS I ONAL TIME 


C DT = 


TIME STEP 


C H = 


NONDIMENS lONAL GRID LENGTH PARAMETER= DX = DY 


C RO = 


NONDIMENS lONAL VORTEX CORE RADIUS 


C XWIDTH= 


NONDIMENS lONAL WIDTH OF COMPUTATIONAL AREA 


C ZTMAX = 


MAXIMUM VORTEX STRENGTH 


C SP = 


NONDIMENS lONAL STRATIFICATION PARAMETER 


C FV = 


NONDIMENS lONAL FROUDE NUMBER 


C BVND = 


NONDIMENSIONAL BRUNT-VALSA FREQUENCY 



C 



C INTEGER 


VARIABLES; 


o 

11 


NUMBER OF NODES IN X DIRECTION 


C N = 


NUMBER OF NODES IN Y DIRECTION 


o 

11 


X NODAL POINT OF VORTEX START POSITION 


o 

11 


Y MODAL POINT OF VORTEX START POSITION. 


n 

It 


X ARRAY DIMENSION 


C NN = 


Y ARRAY DIMENSION 


C NFIG = 


PLOT NUMBER 


C NSTEP = 


ITERATION NUMBER 


C 


\ 



C ARRAY VARIABLES; 



C RHO = 


DENSITY PERTUBATIONS AT EACH NODE DUE TO DYNAMICS 


C RHOI = 


INITAL VALUES OF DENSITY AT EACH NODE DUE TO GRADIENT 


C RHOT = 


SUM OF BOTH DENSITY EFFECTS AT EACH NODE 


C RHNEW = 


NEW VALUE OF RHO AT NEXT TIME STEP 


C ZETA = 


VORTICITY AT EACH NODE 


C ZTNEW = 


NEW VALUE OF VORTICITY AT NEXT TIME STEP 
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C U = X COMPONENT OF VELOCITY AT EACH NODAL POINT 

C V = Y COMPONENT OF VELOCITY AT EACH* NODAL POINT 

C U2 = U VELOCITYS COMPUTED AT BIOT-SAVART NODAL POINTS 

C V2 = V VELOCITYS COMPUTED AT BIOT-SAVART NODAL POINTS 

C X = X COORDINATES AT EACH MODAL POINT 

C Y = Y COORDINATES AT EACH MODAL POINT 

C XU = X COORDINATE AT EACH U 3. S. MODAL POINT 

C YU = Y COORDINATE AT EACH U 3. S. MODAL POINT 

C XV = X COORDINATE AT EACH V B. S. NODAL POINT 

C YV = y' COORDINATE AT EACH V B. S. NODAL POINT 

C 

C 

c 

C SUBROUTINES USED IN THIS PROGRAM: 

C 

C SUBROUTINE COUT - OUTPUTS NODE COORDINATES 

C 

C SUBROUTINE OUTPT - OUTPUTS VARIOUS PROGRAM ARRAYS WHEN CALLED 

C 

C SUBROUTINE CPLOT - CREATES ARRAY OF DENSITY AT EACH POINT 

C AND CALLS DISSPLA CONTOURING ROUTINE TO PLOT. 

C 

C SUBROUTINE ZPLOT 

C 

c . ■ - 

C SUBROUTINE DRPLOT 

C 

c 
c 

C SUBROUTINE 2PL - OUTPUTS VORTICITY CONTOURS AS A PRINTER PLOT. 

C 

C SUBROUTINE VPLOT - CREATES COMPLEX VELOCITY VECTORS AT EACH 

C NODAL POINT AMD THEN CALLS THE FOLLOWING 

C SUBROUTINES TO PLOT THEM ON DISSPLA: 

C A) V2PLOT 



- CREATES ARRAY OF VORTICITY AT EACH POINT AND 
CALLS DISSPLA COUTOURIMG ROUTINE TO PLOT. 

- CREATES ARRAY OF DENSITY PERTUBATIONS AT EACH 
POINT AND CALLS DISSPLA CONTOURING ROUTINE 

TO PLOT. 



32 



C B) AKIN 

C C ) AKAX 

C 

C SUBROUTINE COM - CALLED BY 2PL TO CONTOUR VORTICITY ON PRINTER. 

C 

c 

C FUNCTIONS USED IN THIS PROGRAM: 

C 

C FUNCTION VBS - COMPUTES THE V COMPONENT OF VELOCITY USING 

C THE BIOT-SAVART LAW. 

C 

C FUNCTION UBS - COMPUTES THE U COMPONENT OF VELOCITY USING 
C THE BIOT- SAVART LAW. 

C 

C FUNCTION PNEW - TAKES THE FINITE DIFFERENCE TO STEP THE 
C EQUATIONS IN TIME. 

C 

DIMENSION RHOI(50,100) , U( 50 , 100 ) , V(50 , 100 ) , ZETA( 50 , 100 ) 

* ,2TNEW(50,100) ,RHNEW(50 , 100) ,NPLOT(3) ,X(50) ,Y(100) ,RHO(50 , 100) 

/ 

* ,XU(50) ,YU(100) ,xy(50) ,YV(100) ,U2(50,100) ,V2(50,100) 

CHARACTER'S LABEL 

COMMON WK( 15000) 

COMMON /PN/ DT,I, J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT,FV 
COMMON /PARK/ H;MvM,PI '• ; - 

DATA MAXSTP/40 / 

DATA NPLOT/1 , 10,20,30,40,50,60,70/ 



C 

C ( CALL GRAPHIC OUTPUT DEVICE ) 

CALL COMPRS 
C CALL TEK613 

C 



c 

MM= 50 






MM= 100 



c 



C INPUT PROGRAI^ PARAMETERS 

C 

DT= 2.000E0 
T=0.0 

XWIDTH =5.0 
R0= 0.09E0 
SP= 1.20 
FV=0. 029755 
ZTMIN= l.OE-8 
M= 44 
N= 84 
ISAVE= 1 

C DESIGNATE NODAL POINTS TO PLACE INITAL VORTEX AT (MP,NP) 

MP=6 

NP=66 

C COMPUTE ADDITIONAL PARAMETERS 

C 

C..... DENSITY GRADIENT IN Y DIRECTION INITALLY = DRHO 
DRHOI=-SP*SP*FV*FV 

C NONDIMENS I ONAL B. V. FREQUENCY = BVND 

BVND= SP*SP*FV*FV 

C MAXIMUM VORTICITY AT ANY POINT = ZTMAX 

ZTMAX=-(FV/(R0^R0)) 

C GEOMETRIC CONSTANT PI 

PI= 4.0E0*ATAN(1.0E0) 

C GRID LENGTH PARAMETER = H 

H= XWIDTH/ (M- 4) 

C 

C COMPUTE MISC M AND N PARAMETERS 

MM1= M-1 
MP1= M+1 
MM2= M-2 
NP1=M+1 
NM1=N-1 
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NM2= N-2 
MW = (2’^M) - 7 
NW = NM2 - 1 
C 

C IMITALIZE COORDINATE SYSTEM ORIGINS 

C 

X(l)= -H 
Y(l)= H 
XU(1)= -l.SEO’^H 
YU(1)= H 
XV(1)= -H 
YV(1)= 1.5E0*H 

C 

C COMPUTE COORDINATES 

C ( NOTE: COMPUTATIONAL GRID IS IN 4TH QUADRANT +X -Y ) 
C 

DO 10 1= 2,MP1 

XU(I)= XU(I-1)+H 
XV(I)=XV(I-1)+ H 

10 X(I) = X(I-l) + H 
DO 11 J= 2,NP1 

YU(J) = YU(J-1) - H 
YV(J) = YV(J-l) - H 

11 Y(J) = Y(J-l) - H 
,YHT = ABS(Y(N-2)) 

C OUTPUT INITAL PARAMETERS OF THIS RUN OF PROGRAM 

C 

WRITE(6, 1275) 

WRITE (6, 1500) 

WRITE(6,1535) DT,H 
WRITE(6,1540) XWIDTH,YHT 
WRITE(6,1550) SP,FV 
WRITE(6,1560) M,N 
WRITE (6, 1570) BVND 
C 
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C OUTPUT COORDINATES OF NODES 

WRITE(6,1580) X(MP),Y(NP) 

WRITE(6,1590) 

CALL COUT(X,Y,MM,NN) 

C 

C SPECIFY INITAL VORTICITY AND DENSITY FIELDS 

C 

C ( FOR GRID SYSTEM X,Y COMPUTATIONAL GRID IN THE 4TH QUADRANT ) 
C ( TERMS ARGl - ARG4 ARE IN QUADRANTS 1 - '4 ) 



ARGMN= 30.0 
XCENT= X(MP) 

YCENT= Y(NP) 

DEM= 2.0E0^R0*R0 
ZETA(MP,NP)= ZTMAX 
DO 30 J=1,N 

DO 30 1=1, M 
TRM1=0.0E0 
TRM2=0.0E0 
TRM3=0.0E0 
TRM4=0.0E0 

ARG1= ((X(I)-XCENT)*’^2 + (Y( J)+YCENT)**2)/DEN 
IF (ARGl .GT. ARGMN) GO TO 12 
TRM1= EXP (-ARGl) 

12 ARG2= ( (X(r)+XCENT)’^*2 + (Y( J)+YCENT)’^’^2)/DEN 
IF(ARG2 .GT. ARGMN) GO TO 14 
TRI^2= EXP(-ARG2) 

14 ARG3= ( (X(I)+XCEMT)**2 + (Y( J) -YCENT) **2 ) /DEN 
IF(ARG3 .GT. ARGMN) GO TO 16 
TRM3= EXP(-ARG3) 

16 ARG4= ( (X(I)-XCENT)’^^2 + ( Y( J) -YCENT) ^*2 ) /DEN 
IF(ARG4 .GT. ARGMN) GO TO 13 
TRM4= EXP(-ARG4) 

18 IF((J.EQ.NP) .AND. (I .EQ.MP)) GO TO 20 

ZETA(I , J)=ZTMAX*(-TRM1 +TRM2 -TRM3 +TRM4) 
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IF(ABS(ZETA(I,J)) .LT. ZTMIN) ZETA(I , J)=0 .OEO 
20 RHOI(I,J)= DRHOI*Y(J) 

RHO(I,J)= O.OEO 
30 CONTINUE 

WRITE (6, 1250) 

C CALL OUTPT(ZETA, 'ZTA(I, J) ' ,ra-I,NN,0) 

C CALL OUTPT(RHOI, 'ROI(I, J) ' ,MM,NM,0) 

C 

C ' IMITALIZE TIME STEP COUNTERS 

NFIG=0 

NSTEP=0 

C 

C*x*x*x**** ENTRY POINT OF TIME LOOP 

60 CONTINUE 

C 

C COMPUTE U VELOCITY FOR ALL NODES IN U VELOCITY GRID 

C 

DO 70 1= 1,MP1 
DO 70 J= 1,N 

70 U2(I,J)= UBS(XU(I) ,YU(J) ,MM,NN,X,Y,ZETA) 

C 

C COMPUTE V VELOCITY FOR ALL NODES IN V VELOCITY GRID 

C 

DO 100 1= 1,M 
DO 100 J= 1,NP1 

100 V2(I,J)= VBS(XV(I) ,YV(J) ,MM,NN,X,Y,ZETA) 

C 

C INTERPOLATE VELOCITIES FROM VELOCITY GRIDS TO NODAL 

C POINTS ON COMPUTATIONAL GRID. 

VMIN= l.OE-08 
DO 120 1= 1,M 
DO 120 J= 1,N 

U(I,J)= (U2(I,J) + U2(I+1, J))/2.0E0 
V(I,J)= (V2(I,J)+ V2(I, J+1))/2.0E0 

IF (ABS(U(I,J)) .LT. VMIN) U(I,J)= O.OEO 
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IF (ABS(V(I,J)) .LT. VMIN) V(I,J)= O.OEO 
120 CONTINUE 

C IF (NSTEP .EQ. 0) CALL OUTPT (V, ' V( I , J) ',MM,NN,1) 

C IF (NSTEP .EO. 0) CALL OUTPT(U, ' U(I , J) ',MM,NN,1) 

C 

C 

C 

C INCREASE TINE COUNTER 3Y ONE. 

C CHECK ITERATION COUNTER AND TERMINATE IF NSTEP IS MORE THAN MAX 
C COMPUTE VORTICITY AND DENSITY AT ALL INTERIOR POINTS 

C ( STORED AS ELEMENTS OF ZTNEW AND RHNEW ) 

C 

c 

NSTEP=NSTEP+1 

IF (NSTEP .GT. MAXSTP) GO TO 1700 
T= T+DT 

DO 600 1=2, MMl 

DO 600 J=2,NM1 
VPT= V(I,J) 

UF= (U(I+1,J)+U(I, J))/2.0E0 
UFA=ABS(UF) 

UB= (U(I-1, J)+U(I,J))/2.0E0 
UBA=ABS(UB) 

VF= (V(I, J-1)+V(I, J))/2.0E0 
VFA=ABS(VF) 

VB= (V(I,J+1)+V(I,J))/2.0E0 
VBA=ABS(VB) 

ZTNEW ( I , J ) = PNEW ( ZETA , RHO , MM , NN , - 1 . OEO , 0 . OEO , 0 . OEO ) 



600 RHNEW(I,J)= PNEW(RHO,RHO,MM,MN,O.OEO,O.OEO,BVMD) 

C 

C 

C APPLY BOUNDRY CONDITIONS TO_ RHO AND ZETA 

C 

C ASSIGN NEW VALUES OF RHO AMD ZETA TO RHO AND ZETA ARRAYS . . 

C 
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c 



DO 700 I = 2,MM1 
DO 700 J = 2,NH1 

2ETA(I,J)= ZTMEW(I,J) 

700 RHO(I, J)=RHMEW(I , J) 

C CALCULATE RHO AND ZETA VALUES AT I=M-2 AMD J=N-2 BOUNDRIES ... 

C 

DO 300 1= 2,HM2 

DO 300 J= 2, MHZ 

ZETA(MH2,J)= (ZETA(MM1,J)+ ZETA(M-3, J))/2.0E0 
ZETA(I,NH2)= (ZETA(I,NM1)+ ZETA( I ,N-3 ) ) /2 . OEO 
RHO(MM2,J)= (RH0(MH1,J)+ RHO (M-3 , J) ) /2 . OEO 
RHO(I,NH2)= (RH0(I,NH1)+ RHO ( I ,N-3 ) )/2 . OEO 

IF(ABS(ZETA(I,J)) .LT. ZTMIN) ZETA(I , J)=0 .OEO 
IF(ABS(RHO(I, J) ) .LT. ZTMIN) RHO(I,J)= O.OEO 
BOO CONTINUE 
C 

C 

C PLOT FLOW PATTERN WHEN NSTEP = VALUE SPECIFIED IN NPLOT 
C 

C THEN RESTART ALGORITHM AT LINE 60 FOR A NEW TIME STEP 

C. 

DO 1100 K= 1,3 

IF (NSTEP-NPLOT(K) ) 1100, 1200, 1100 

1100 CONTINUE 
GO TO 60 
1200 NFIG=NFIG+1 
TM= T’^FV 

IF (MSTE? .EQ. HAXSTP ) CALL OUTPT(ZETA, ' ZTA(I ,J) ', MM, NN, NSTEP) 
IF (NSTEP .EQ. MAXSTP) CALL OUTPT (RHO, • ' RHO ( I , J) , MM ,NN ,NSTEP ) 
WRITE (6, 1300 )NSTEP,T,TM 

CALL CPLOT ( RHO , RHOI , MM , NM , NSTEP , TM , MW , MW ) 

CALL VPLOT ( U , V , X , Y , MM , NN , TM ) 

CALL DRPLOT ( RHO , MM , NM , NSTEP , TM , MW , NW ) 

CALL ZPL(ZETA,MM,NN,XWIDTH,YHT, ITER) 
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CALL 2PL0T(ZETA,MM,NN,TM,MW,NW) 

1250 FORMAT(// 10X,'INITAL ZETA AND RHO VALUES : ') 

1275 FORMAT( ' 1 ' ) 

1300 FORI-IAT(/10X, ' ITERATION MUM : ',16 , 4X , ' TIME= ' , F9 . 3 , 4X , ' T’^ =’,F10.6) 
1400 WRITE (6, 1500) 

1500 FORMAT ( ///) 

1530 FORMATC ',5X,'INITAL PROGRAM PARAMETERS:' ) 

1535 FORMAT( 'O' ,5X, 'TIME STEP= ' , ?3 . 4 , 5X , ' GRID SPACING H = ',”8.4) 
1540 FORMAT( 'O' ,5X, 'NON DIMENSIONAL CELL DIMENSIONS : XDIRECTION =', 

* F8.4,5X, 'YDIRECTION =',F8.4) 

1550 FORMAT('0' ,5X, 'STRATIFICATION PARAMETER =',F8.4,5X, 

* 'FROUDE NUMBER =',F8.4) 

1560 FORMAT( 'O' ,5X, 'NUMBER OF NODES IN X DIRECTION =',I4,5X, 

* ' IN Y DIRECTION =' ,14) 

1570 FORMAT( 'O' ,5X, 'NONDIMENSIONAL 3. V. FREQUENCY =',F15.9) 

1580 FORMAT( 'O' ,5X, 'X COORDINATE OF VORTEX START POSITION =',F8.4,5X, 
'Y COORDINATE =',F8.4) 

1590 FORMAT( 'O' ,5X, 'OUTPUT COORDINATES OF ALL OTHER NODES:') 

C 

C LOOP BACK TO STARTING POINT OF ITERATION 

GO TO 60 

1700 CALL DONEPL 
C1700 CONTINUE 
C 

C SAVE ARRAYS AND TIME DATA TO RESTART PROGRAM LATER 

C 

IF (ISAVE .EQ. 0 ) GO TO 1900 
WRITE(3,*) SP,FV,DT,T 
WRITE ( 8 , ^ ) M , N , MP , NP 
WRITE (8,^) BVND,H,ZTMAX 
DO 1800 1= 1,M 

DO 1800 J= 1,N 

1800 WRITE(8,*) ZETA(I,J), RH0(I,J) 

1900 STOP 
END 
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c 




REAL FUNCTION PNEW (P , Q, MM , NN , A , B , C) 

C 

C THIS FUNCTION INTEGRATES THIS EQUATION WITH RESPECT TO TIME 
C $$$ DP/DT= (-D(UP)/DX -D(VP)/DY + A*DQ/DX + B’^DEL(P) +C*V ) $$$ 

C UTILIZING AN UPWIND DIFFERENCING SCHEME ' 

C 

DIMENSION P(HM,ra),Q(MM,NN) 

COMl^ON /PN/ DT,I, J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT,FV 
COra-ION /PARM/ H,M,N,PI 

Pl= (UF-UFA)*P(I+1,J) +(UF+UFA-UB+UBA)’^P(I, J) - (UB+UBA)’^P(I-1 , J) 
P2= (VF-VFA)*P(I, J-1) +(VF+VFA-VB+VBA)*P(I, J) - (VB+VBA)*P ( I , J+1 ) 
P3= Q(I+1,J) -Q(I-1,J) . 

P4f P(I+1,J) +P(I-1,J)+P(I,J-1)+P(I,J+1)-4.0E0’^P(I, J) 

P5= 2.0E0*H*VPT 

TEMP1=(DT/(2.0E0*H))*((-P1-P2) +A*P3+ (2 . OEO/H) *B*P4+C*P5 ) 
PNEW=P(I, J)+TEMP1 

• RETURN ... . • ■ 

END 

C 

C==================================================================== 

c 

REAL FUNCTION ' UBS ( XU , YU , MM , NN , X , Y , ZETA ) 

C 

C THIS SUBROUTINE FINDS FOUNDRY VALUES OF VELOCITY U 
C UTILIZING THE BIOT-SAVART LAW 

C 

C (NOTE: TERMS IN EQUATION ARE EVALUATED CONSECUTIVELY IN EACH 
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c 



QUADRANT FROM ONE TO FOUR ) 



C 

DIMENSION X(MM) , Y(NN) , ZETA(MM,NN) 

COMMON /FARM/ K,M,N,PI 
C 

UHIN= l.OE-08 
Ul= O.OEO 
U2= O.OEO 
U3= O.OEO 
U4= O.OEO 
C 

SIGNP=1 .OEO 
SIGNM=-1 .OEO 
DA= H*H 
C 

C - COMPUTE U COMPONENT OF VELOCITY USING THE B-S LAW 

c - ■ 

DO 20 1= 1,M 

DO 20 J= 1,N 
C -- QUADRANT ONE -- 

RSQ1=((XU-X(I))**2 + (YU+Y(J))**2) 

Ul= U1+ ((-Y(J)-YU)/(2.0E0*PI*RSQ1))*SIGNM*2ETA(I,J)*DA 
C -- QUADRANT TWO -- 

RSQ2= ((XU+X(I))**2 +(YU+Y(J))*’^2) 

U2= U2+ ((-Y(J)-YU)/(2.0E0^PI*RSQ2) )’^SIGNP*2ETA(I, J)*DA 
C -- QUADRANT THREE -- 

RSQ3= ((XU+X(I) )*’^2 +(YU-Y(J) )**2) 

U3= U3+((Y(J)-YU)/(2,0E0’^PI’^RSQ3) )’^SIGNM^ZETA(I, J)*DA 
C -- QUADRANT FOUR -- 

RSQ4= ((XU-X(I))’^^2 + (YU-Y(J) )’^*2) 

20 U4 = U4 + ((Y(J)-YU)/(2.0E0*PI*RSQ4))*SIGMP*ZETA(I,J)’^DA 
UTEMP= U1+U2 +U3+ U4 
IF (ABS(UTEHP) .LT.UHIN) UTEMP=0.0E0 
UBS= UTEMP 
RETURN 
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END 



C 

C 

c================================================================== 

c 

REAL FUNCTION VBS (XV , YV ,MM ,NN , X, Y , ZETA) 

C 

C. THIS SUBROUTINE FINDS BOUMDRY VALUES OF VELOCITY V 
C UTILIZING THE BIOT-SAVART LAW 

C 

C (NOTE : TERNS IN THE EQUATION ARE EVALUATED CONSECUTIVELY IN EACH 
C QUADRANT FROM ONE TO FOUR) 

C 

DIMENSION X(MM) ,Y(NN) ,ZETA(MN,NN) 

COMMON/PARM/ H,M,N,PI 
C 

VMIN=: l.OE-08 
Vl= O.OEO 
V2= O.OEO 
V3= O.OEO 
. V4= O.-OEO 

C ' . ■ 

SIGNP=1.0E0 
SIGNN=-1 .OEO 

.DA= H*H • • 

C 

c COMPUTE V COMPONENT OF VELOCITY USING THE B-S LAW 

C 

DO 20 1= 1,M 

DO 20 J= 1,N 
C -- QUADRANT ONE -- 

RSQ1=((XV-X(I))**2 + (YV+Y(J))**2) 

Vl= V1+ ((XV-X(I))/(2.E0*PI^RSQ1))*SIGNT'I*ZETA(I,J)*DA 
C -- QUADRANT TWO -- 

RSQ2= ( (XV+X(I) )**2 +(YV+Y( J) )*’^2) 
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V2= V2+ ((XV+X(I))/(2.0E0’^PI*RSQ2))*SIGNP’^ZETA(I,J)’^DA 
C. -- QUADRANT THREE -- 

RSQ3= ((XV+X(I))’^’^2 +(YV-Y(J))’^*2) 

V3= V3+ { (XV+X ( I ) ) / ( 2 . 0E0*PI’^RSO3 ) ) *SIGMM’^ZETA( I , J ) ’^DA 
C -- QUADRANT FOUR -- 

RSQ4= ( (XV-X(I) )’^’^2 + (YV-Y(J))^’^2) 

20 V4 = V4 + ((XV-X(I))/(2.0E0*PI*RSQ4))*SIGNP’^ZETA(I,J)*DA 
VTEMP = V1+V2 +V3+ V4 
IF(ABS(VTSMP) .LT.VNIN) VTENP=0.0E0 
VBS= VTENP 
RETURN 
END 
C 

c 

C 

c 

SUBROUTINE COUT(X,Y,MN,NM) 

C 

c 

C , 

REAL X(MM) ,Y(NN) ,H,PI 
COMMON/PARM/H , M , N , PI 
WRITE(6,50) 

WRITE (6,*) ' ARRAY X' 

■ •• WRITE(6,*) 

WRITE(6,*) (X(I), 1=1, M) 

WRITE(6,*) 

WRITE(6,*) 

WRITE (6,*) ' ARRAY Y 

WRITE (6,*) 

WRITE(6,*) (Y(J), J=1,N) 

WRITE(6,^) 

WRITE (6,*) 

50 FORMAT (///) 

RETURN 
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END 



C===================================================================: 

SUBROUTINE OUTPT ( A , LBL , MM , NN , KAL ) 

C 

C THIS SUBROUTINE PRINTS OUT AN (NKM) ARRAY WITH A LABEL 
C IN A RECTANGULAR ( 2 , M- 2 )X( 2 ,N-2 ) GRID 

C 

DIMENSION A(MM,MN) 

COMMON/PARM/ K , M / N , P I 
CHARACTER’^8 LBL ' 

NM2= N-2 
MM2= M-2 
WRITE (6, 50) 

WRITE (6,*) 

WRITE (6,*) ' OUTPUT OF ARRAY: ',LBL,' ITER NUM KAL 

WRITE (6,*) 

DO 10 J= 2,NM2 
WRITE (6,=^) 

WRITE(6,*) (A(I,J), 1=2, MM2) 

10 CONTINUE 
50 FORMAT( '1' ) 

RETURN 

END 

C 

C- - - 

c 

SUBROUTINE CPLOT ( RHO , RHOI , MM , NN , NSTEP , TM , MW , NW) 

COMMON WK( 15000) 

COMMON /FARM/ H,M,N,PI 
REAL W(81,81) 

REAL RHO(MM,NN) ,RHOI(MM,NM) , H , PI , RHOT( 50 , 100) 

DATA ISCALE/4HSCAL/ 

C 

C... CREATE ARRAY OF DENSITY CHANGES AT EACH NODE (DRHO) 

C 
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MM2= M-2 
MM2= N-2 
DO 10 1= 2,HM2 

DO 10 J= 2,NH2 
JJ= (NM2+1)-J 
11= I+(M-5) 

RHOT ( I , J ) =RHO ( I , J ) +RHO I ( I , J ) 
W(II,JJ)= R>:OT(I,J) 

10 CONTINUE 

DO 20 1= 3, MM2 
11= (M-l)-I 

DO 20 J= 2,NM2 
JJ= (NM2+1) -J 



20 W(II,JJ)= RHOT(I,J) 

C 

C CALL OUTPT ( RHOT , ' RHT ( I , J ) J , MM , NN , NSTEP ) 

C CALL SUBROUTINE CONTR TO PLOT THESE PERTU3ATI0NS 

C 



CALL RESET (3HALL) 

CALL PAGE(8.7,11.2) 

CALL PHYSOR(2.0,3.00) 

CALL AREA2D (5. 0,6.0) 

CALL HEIGHT (.2 ) 

CALL INTAXS 

‘CALL'XNAME ('X/BO',4) 

CALL YNAME ( 'Y/BO' ,4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF ( - 5 . 0 , 1 . , 5 . , 10 . , 1 . , 0 . ) 

CALL MIXALF( ' INSTRU' ) 

CALL MESSAG('T(EH.5)^(EXHX) =$',100,1.5,6.5) 
CALL REALNO(TM, 106 , 'ABUT' , 'ABUT' ) 

CALL FRAME 

CALL BCOMON (15000) 

CALL CONMIN(2.0) 
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CALL CONANG (20.) 

CALL CONMAK ( W , MW , MW , ISCALE ) 

CALL HEIGHT (.08) 

CALL CONLIN ( 0 , 5HSOLID , ' LABELS ' 2 , 5 ) 

CALL CONLIN ( 1 , 4HDASH , 3HNOLABELS , 1 , 3 ) 

CALL RASPLN(0.25) 

CALL COMTUR ( 2 , 6HLABELS , 4HDRAW) 

CALL HEIGHT (.2) 

CALL COMPLX 

CALL RESET ( 'HEIGHT' ) 

CALL RESET ('COMPLX') 

CALL ENDPL (0) 

RETURN 
END ' 

C 

C - 

C 

SUBROUTINE ZPLOT ( ZETA , MM , NN , TM , MW , NW ) 

REAL ZETA(MM,NN) ,H,PI 
COMMON WK( 15000) 

COMMON /PARM/ H,M,N,PI . , 

DIMENSION W(81,81) 

DATA ISCALE/4HSCAL/ 

C 

C... CREATE ARRAY FOR PLOTTING OF ZETA IN ACTUAL COMPUTATIONAL AREA.. 
C 

m2= M-2 
NM2=N-2 

DO 10 1= 2, MM2 
11= I+(M-5) 

DO 10 J= 2,NM2 
JJ= (NM2+1)-J 
10 W(II,JJ)=ZETA(I,J) 

DO 20 1= 3, MM2 
11= (M-l)-I 
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c 

c 

c 

c 



DO 20 J= 2,NM2 
JJ= (N112+1) - J 
W(II,JJ) = -ZETA(I,J) 

CALL SUBROUTINE COMTR TO PLOT VORTICITY 



CALL RESET (SHALL) 

CALL PAGE (8. 7, 11. 2) 

CALL PKYSOR(2.0,3.00) 

CALL AREA2D (5. 0,6.0) 

CALL HEIGHT (.2 ) 

CALL INTAXS 
CALL XMAME ( ' X/BO ' ,4) 

CALL YNAllE ( ‘Y/BO ' ,4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF ( -5 . 0 , 1 . , 5 . , 10 . , 1 . , 0 . ) 

CALL MIXALF( ‘INSTRU' ) 

CALL MESSAG( ’T(EH.5)’^(EXHX) =$',100,1.5,6,5) 
CALL REALNO(TM,106, 'ABUT' , 'ABUT' ) 

CALL FRAME 

CALL BCOMOM (15000) 

CALL CONAMG (60.) 

CALL CONMIN(2.0) 

CALL CONMAK (W,MW,NW, ISCALE) 

CALL HEIGHT (.03) 

CALL CONLIN (0 , 5HSOLID , ' LABELS ' , 2 , 5 ) 

CALL CONLIN ( 1 , 4HD ASH ,' LABELS ', 1 , 3 ) ' 

CALL RASPLN(0.25) 

CALL CONTUR ( 2 , 6HLABELS , 4HDRAW) 

CALL HEIGHT ( .2) 

CALL COMPLX 

CALL RESET ( 'HEIGHT' ) 

CALL RESET ( ' COMPLX' ) 
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CALL ENDPL (0) 

RETURN 

END 

C 

C 

C-- - 

c 

SUBROUTINE DRPLOT ( RKO , MM , MM , MSTE? , TM , MW , MW ) 

COMMON WK( 15000) 

COmON /FARM/ H,M,N,PI 
REAL W(81,81) 

REAL RHO(MM,NN) ,H,PI 
DATA ISCALE/4HSCAL/ 

C 

C... CREATE ARRAY OF DENSITY CHANGES AT EACH NODE (DRHO) 
C 

MM2= M-2 
NM2= N-2 
DO 10 1= 2, MM2 

DO 10 J= 2,NM2 
JJ= (NM2+1)-J 
11= I+(M-5) 

W(II,JJ)= RHO(I,J) 

10 ' CONTINUE 

DO 20 I=, 3, MM2 
■ 11= (M-l)-I 

DO 20 J= 2, MM2 
JJ= (NM2+1) -J 



20 W(II,JJ)= RHO(I,J) 

C 

C CALL SUBROUTINE CONTR TO PLOT THESE PERTUBATIONS 

C 



CALL RESET (SHALL) 
CALL PAGE(8.7,11,2) 
CALL PHYSOR(2.0,3.00) 
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CALL AREA2D (5. 0,6.0) 

■ CALL HEIGHT( .2 ) 

CALL IMTAXS 

CALL XNAME ('X/BO' ,4) 

CALL YMAME ( ' Y/BO ' ,4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL .GRAF (-5.0,1.,5.,10.,1.,0.) 

CALL MIXALF( ' INSTRU' ) 

CALL MESSAG('T(EH.5)*(EXHX) =$',100,1.5,6.5) 

CALL REALNO (TM , 106 , ' ABUT ' , ' ABUT ' ) 

CALL FRAME 

CALL BCOMOM (15000) • 

CALL CONMIN(2.0) 

CALL CONAMG (20.) 

CALL CONMAK (W,MW,^^■/, ISCALE) 

CALL HEIGHT (.08) 

CALL CONLIN (0 , 5HSOLID , ' LABELS ' , 2 , 5 ) 

CALL CONLIN ( 1 , 4HDASH , 8HNOLABELS , 1 , 3 ) 

CALL RASPLN(0.25) 

CALL CONTUR ( 2 , 6HLABELS , 4HDRAW) 

CALL HEIGHT (.2) 

CALL COMP LX 

CALL RESET ('HEIGHT') 

CALL RESET ( 'COMPLX' ) 

CALL ENDPL (0) 

RETURN 

END 

C 

C====================================================================== 

c 

SUBROUTINE CON ( 22 , XRANGE , YRAMGE , NX , H , MM , NN , ZMAX , ZMIN ) 

C 

C THIS SUBROUTINE GENERATES A CONTOUR PLOT OF ZETA IN A RECTANGULAR 

C DOMAIN OF SIZE: YRANGE * XRANGE ON THE PRINTER. 
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c 



CHARACTER’^ 1 SYMBOL (8) , GRAPH (i 00) 

REAL LX,LY 

DIMENSION Z2(MM,NN) , VALUE (7) 

DATA SYl-IBOL / ' X ' , ' A ' , ' - ' , ' ' , ' ' , ' + 7 

VALUE (1)= ZMIN 

VALUE(2)= ZMIN*. 7 5 

VALUE (3 )= ZMIN*. 5 

VALUE(4)= 0.0 

VALUE (6 )= ZMAX*.75 

VALUE (7 )= ZMAX 

VALUE(5)=ZMAX*.45 

NY= 0 . 8*NX* YRANGE / XRANGE 

DELX= XRANGE/NX 

DELY= YRANGE/NY 

WRITE (6, 50) 

WRITE ( 6 , * ) ' NX= ' , NX , ' NY= ' , NY , ' XRANGE= ' , XRANGE , 

* ' YRANGE= ' , YRAI'IGE 

WRITE (6,*) 

DO 11 1= 1,NX 
11 GRAPH(I)= 'T' 

WRITE(6,6) (GRAPH(I), 1= 1,NX) 

DO 5 JSYMBL = 1,NY 
Y= (JSYMBL-0.5)*DELY 
J= 1+Y/H 
LY= Y-(J-1)*H 
HMLY= H-LY 

DO 4 ISYl^BL = 1,NX 

X= (ISYM3L-0.5)*DELX 

I=l+X/H 

LX= X-(I-1)*H 

HMLX= H- LX 

Al= HMLX*HMLY 

A2= HMLX*LY 

A3= LX*LY 
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A4= LX’^HMLY 

Z2C= (A1*Z2(I,J) +A2’^Z2(I, J+1) +A3*Z2(I+1, J+1) 

* +A4=*'Z2(I+1, J))/H*’^2 
C 

C... DETERMINE THE VALUE OF MKAMGE BASED OM Z2C 

C 

C PRINT OUT CHARACTERS IN ARRAY GRAPH BY ROW 

C 

DO 2 K= 1,7 

IF (Z2C-VALUE(K)) 1, 1, 2 

1 NRANGE = K 

GO TO 3 

2 CONTINUE 
NRANGE = 8 

3 GRAPH ( ISYMBL)= SYMBOL (NRANGE ) 

4 CONTINUE 

5 WRITE(6,6) (GRAPH(I), 1=1, NX) 

6 FORMAT (' ',10X,100A1 ) 

50 FORMAT ( ' 1 ' ) 

DO 60 1= 1, NX 
60 GRAPH(I)= 'B' 

. WRITE(6,6) (GRAPH(I) ,1 = 1, NX) . 

WRITE(6,*) 

WRITE (6,*) 'MAXIMUM VALUE OF NEGATINE VORTICITY = SYMBOL (X) = 

* ,ZMIN . . - 

WRITE(6,*) ' 75% OF MAXIMUN NEGATIVE VORTICITY = SYMBOL (A) ' 
WRITE (6,^) ' 50% OF MAXIMUM NEGATIVE VORTICITY = SYMBOL (-) ' 
WRITE (6,*) 

WRITE (6,*) ' MAXIMUM VALUE OF POSITIVE VORTICITY = SYMBOL (1) = 
, ZMAX 

WRITE(6,^) ' 75% OF MAXIMUM POSITIVE VORTICITY = SYMBOL (+) ' 

RETURN 

END 

C 

c 
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SUBROUTINE ZPL (ZETA , MM , NN , XWIDTH , YHT , ITER ) 

COMMON /FARM/ H,M,N,PI 
DIMENSION Z2(50,100) ,ZETA(MM,MN) 

C 

C... CREATE ARRAY FOR PLOTTING OF ZETA IN ACTUAL COMPUTATIONAL AREA.. 
C 

MM3= M-3 . 

MM3=N-3 
ZMIN= 0.0 
ZMAX= 0.0 
DO 10 1= 1,MM3 

DO 10 J= -1 ,NM3 
Z2(I,J)= ZETA(I+1,J+1) 

C.. FIND THE MAXIMUM POSITIVE VORTICITY IN ARRAY Z2(MM,NN) 

C 

IF(Z2(I,J) .GT. ZMAX) ZMAX=Z2(I,J) 

C.. FIND THE MAXIMUM NEGATIVE VORTICITY IN ARRAY Z2(MM,NN) 

C 

IF( Z2(I,J) .LT. ZMIN) ZMIN= Z2(I,J) 

10 ■ CONTINUE 

C 

C ‘ • 

IF (ZMAX .LT. l.OE-6 ) ZMAX=1.0 

C CALL SUBROUTINE CON TO PLOT VORTICITY 

c ■ 

CALL CON ( Z2 , XWIDTH , YHT , 1 0 0 , H , MM , NN , ZMAX , ZMIN ) 

RETURN 

END 

C- 

SUBROUTINE VPLOT ( U , V , X , Y , MM , NN , TM ) 

COMMON /FARM/ H,M,N,PI 
COMPLEX ZV( 50, 100) , CV(50,100) 

DIMENSION U(MM,NM), V(MM,NN), X(MM) , Y(NN) 

MM2= M-2 
NM2= N-2 
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MM3 = M - 3 
MM3 = N - 3 
C 

C FORM COMPLEX U, V, AND Z ARRAYS FOR QUADRANTS 3 AND 4 

C 

DO 10 1= 2, MM2 

DO 10 J= 2, MM2 
11= I-l 

JJ= (NH2+1) -J 
ZV(II,JJ)= CMPLX(X(I) ,Y(J) ) 

10 CV(II,JJ)= CMPLX(U(I,J),V(I,J)) 

C 

C 

C PLOT VECTORS USING SUBROUTINE V2PLOT 

CALL V2PLOT(ZV,CV,MM3,NM3,TM) 

C 



RETURN 

END 



C 

C 

C • 

SUBROUTINE V2PLOT (Z , C , M ,N ,TM) 



C 

C 

C 

C 

C 

C 



A X ^ A X X X X A A X X A 7C X A X A A X X X X X X X A A A X 7^ X 

* SUBROUTINE V2PLOT TO PLOT MXN COMPLX VELOCIT FIELD C(M,N) ^ 

* OVER A COMPLX DOMAIN Z(M,N) * 

t: a 






COMPLEX Z(50,100) ,C(50,100) ,CMAX 
REAL X(IOO) ,Y(100) ,YX(100) ,YN( 100) ,XX(100) , 
* XN(IOO) ,X2(100) ,Y2(100) 

INTEGER PMAX,PMIN 
EXTERNAL AMAX 
EXTERNAL AMIN 



C 
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DO 20 J =1,N 



DO 10 I =1,M 

X(I) =ABS(REAL(C(I, J))) 
Y(I) =ABS(AIMAG(C(I,J))) 
10 CONTINUE 

CALL AMAX(Y,M,YMAX,PMAX) 
CALL AMAX(X,M,XI'IAX,?MAX) 
YX(J) = YMAX 
XX(J) = XMAX 
CALL AMIN(Y,N,YNIN,PMIN) 
CALL AMIN(X,M,XMIN,PMIN) 
YN(J) =YNIN 
XN(J) =XMIN 
20 _ CONTINUE 

• CALL AMAX(YX,N,YMAX,PMAX) 

C CALL AMIN(YN,N,YNIN,PMIN) 

CALL AMAX(XX,N,XMAX,PMAX> 

C CALL AMIN(XN,N,XHIN,PHIN) 

C 

CMAX = CMPLX( XMAX, YMAX) 

UMAX =SQRT( (XMAX) ^^^^2+ (YMAX) **2) 
DO 24 I =1,M 
DO 22 J=1,M 

C(I,J) = C(I,J)/UMAX 
22 CONTINUE 

24 CONTINUE 

C 

DO 40 J =1,N 

DO 30 I =1,M 

X(I) =REAL(2(I,J)) 

Y(I) =AIMAG(2(I,J)) 

30 CONTINUE 

CALL AMAX ( Y , M , YMAX , PMAX ) 
CALL AMAX (X,M, XMAX, PMAX) 
YX(J) = YMAX 
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XX(J) = XMAX 

CALL AMIN(Y,M,YMIN,PMIN) 

CALL AMIN(X,M,XMIN,PMIN) 

YH(J) =Y!-IIM 
XN(J) =XHIN 
40 CONTINUE 

CALL AMAX ( YX , N , YMAX , PMAX ) 

CALL AH I N ( YN , N , YH I N , PHI N ) 

CALL AHAX(XX,N,XHAX,PHAX) 

CALL AMIN(XN,N,XHIN,PMIN) 

C 

CALL RESET (SHALL) 

CALL PAGE (8. 7, 11. 2) 

CALL PHYSOR(2.0,3.0) 

CALL AREA2D(5. ,6.) 

CALL HEIGHT (.2) 

CALL INTAXS 

CALL XNAHE( 'X/BO ' ,4) 

CALL YNAHE( 'Y/BO' ,4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF(0 . 0 , 0 . 5, 5. 0,-10. 0 , 1 . 0 , 0 . 0 ) 

CALL HIXALF( ' INSTRU' ) 

CALL HESSAG( 'T(EH.5)*(E:^.X) =$’,100,1.5,6.5) 
CALL REALNO(TM, 106, 'ABUT' , 'ABUT' ) 

C 

DO 60 J=1,N 

DO 50 1=1, M 

X(I) = REAL(Z(I,J)) 

X2(I) = X(I)+REAL(C(I, J)) 

Y(I)= AIHAG(Z(I,J)) 

Y2(I) = Y(I)+AIHAG(C(I,J)) 

50 CONTINUE 

CALL HARKER (3) 

CALL SCLPIC( .001) 
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54 

60 



C 

C 

C 

C 

C 

C 

c 

c 

c 

c 



c 



5000 

5100 



CALL CURVE (X,Y,M,-1) 

DO 54 I =1,M 
. X11=X(I) • 

X12=X2(I) 

Y11=Y(I) 

Y12=Y2(I) 

IF(X12.EO.O. .AND.Y12.EQ.0) GO TO 54 
CALL RLVEC (Xll ,Y11 ,X12 ,Y12 , 1121) 
CONTINUE 
CONTINUE 

. CALL RESET ( 'HEIGHT' ) 

CALL ENDPL(O) 

RETURN 

END 



SUBROUTINE AMAX(Y,N, YMAX,PHAX) 

7*: X 7^ X x X 7^ 7^ X 7C 7^ X 7^ X 7^ 7^ 7^ 7^ X X X 

* SUBROUTINE TO COMPUTE THE LARGE ELELEMENT 

* IN A GIVEN ROW OR COLUMN IN ARRAY ' A ' . 

7^: 7*: T^r 7^: 7^ 7^: 7^: 7^: X 7^ 7k 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^: 7^ 7^ 7^ 7^ 7^ 7^: 7^ 7^ 7^: 7^: 7^: 7^: 7^ 7^: X 7^ 7^: 7^: 

*** VARIABLE DECLARATION 
REAL AMAX ‘ . 

INTEGER PMAX 
DIMENSION Y(N) 

YMAX=Y(1) 

DO 5100 1=2, M 

IF(Y(I) .LE.YMAX) GO TO 5000 
YMAX=Y(I) 

PMAX=I 
CONTINUE • 

CONTINUE 
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c 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



6000 

6100 



MAXROW=POS 

RETURN 

END 



SUBROUTINE ANIN ( Y , N , YMIN , PMIN ) 



X X X 7C X X 7C 7? X X 7? t: X X 7^ X X X X X X t: X 7T ^ 



* SUBROUTINE TO COMPUTE THE SMALL ELEMENT 

* IN A GIVEN ROW OR COLUMN IN ARRAY ' A ' . 

7^ 7^ 7^ 7^ 7^ 7^ 7^ ^ 7^ ^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ 7^ T: 7^ 7^ 7^ X X 7^ 7^ 7^ 



*** VARIABLE DECLARATION 
REAL AMIN 

INTEGER PMIN 
DIMENSION Y(N) 

YMIN=Y(1) 

DO 6100 1=2, N 

IF(Y(I) .GE.YMIN) GO TO 6000 
YMIN=Y(I) 

PMIN=I 

CONTINUE 

CONTINUE 

RETURN 

END 
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Figure B.la Velociiy Field: SP = 0.00, T* = 0.05952. 
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Figure B.lb Vorticity Field; SP = 0.00. 
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Figure B.2a Velocity Field: SP = 0.00. 
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Figure B.2b Vorticity Field; SP = 0.00. 
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Figure B.3a Velocity Field: SP = 0.00. 
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Figure B.3b Vorticity Field: SP = 0.00. 
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Figure B.4a Velocity Field: SP = 0.00. 
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Figure B.4b Vorticity Field; SP 
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Figure B.5a Velocity Field; SP = 0.00. 
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Figure B.5b Vorticity Field: SP = 0.00. 



68 



2=97600 





Figure B.6a Velocity Field: SP = 0.00. 
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Figure B.6b Vorticity Field; SP = 0.00. 
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Figure B.7a Velocity Field: SP 
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Figure B.7b Vorticity Field: SP = 0:00. 
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Figure B.8a Velocity Field; SP = 0.00. 
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Figure B.8b Vorticity Field; SP = 0.00. 
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Figure B.9a Velocity Field: SP = 0.00. 
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Figure B.9b Vorticity Field; SP = 0.00. 
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Figure B.lOa Velocity Field: SP = 0.00. 
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Figure B.lOb Vorticity Field: SP = 0.00. 
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Figure B.lla Velocity Field: SP = 0.00. 
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Figure B.llb Vorticity Field: SP = 0.00. 
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Figure B.12a Velocity Field: SP = 0.00. 
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Figure B.l2b Vorticity Field; SP = 0.00. 
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Figure B.13a Velocity Field: SP = 0.00. 
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Figure B.13b Vorticity Field: SP 
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Figure B.14 Experimental and Numerical Results for Various Initial Core Radii. 
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Fieure B.15a Constant Densitv Contours: SP = 1.00. 
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Figure B.15b Velocity Field: SP 
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Figure B.15c Density Perturbation Contours: SP = 
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Figure B.l5d Vorticity Field; SP = 1.00. 
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Figure B.16a Constant Density Contours; SP 
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Figure B.16b Velocity Field: SP = 1.00. 
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Figure B.16c Density Perturbation Contours: SP = 1.00. 
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Figure B. 16 d Vorticity Field: SP = l.OO. 
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Figure B.17a Constant Density Contours: SP = 1.00. 
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Figure B.17b Velocity Field: SP = 1.00. 
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Figure B.17c Density Perturbation Contours: SP = 1.00. 
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Figure B.17d Vorticity Field: SP = 1.00. 
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Figure B.lSa Constant Density Contours: SP = 1.00. 
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Figure B.lSb Velocity Field: SP = 1. 00. 
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Fiaure B.18d Vorticitv Field: SP = 1.00. 
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Figure B.l9a Constant Density Contours; SP = 1.00. 
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Figure B.I9b Velocity Field: SP = 1.00. 
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Figure B.19c Density Perturbation Contours: SP = 1.00. 
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Figure B.19d Vorticity Fieid: SP = 1.00. 
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Figure B.20b Velocity Field; SP = 1.00. 



107 



Y/BO 



1 .71900 



T 







- 5 - 4 - 3 - 2-10 1 2 3 4 



5 



X/BO 



Figure B.20c Density Perturbation Contours; SP = 1.00. 
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Figure B.20d Vorticity Field; SP = 1.00. 
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Figure B.21a Constant Density Contours; SP = 1.00. 
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Figure B.21b 



Velocity Field: SP = 1.00. 
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Figure B.21c Density Perturbation Contours: SP = 1.00. 
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Figure B.21d Vorticity Field: SP == 1.00. 
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Figure B.22a Constant Density Contours: SP = 1.00. 
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Figure B.22b Velocity Field: SP = 1.00. 
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Figure B.22c Density Perturbation Contours: SP = l.OO. 
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Figure B.23a Constant Density Contours: SP = 1.00. 
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Figure B.23b Velocity Field: SP = 1.00. 
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Figure B.23c Density Perturbation Contours: SP = 1.00. 
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Figure B.23d Vorticity Field: SP = 1.00. 
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Figure B.24a Constant Density Contours: SP = 1.00. 
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Figure B.24b Velocity Field: SP = 1.00. 
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Figure B.24c Density Perturbation Contours: SP = 1.00. 
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Figure B.24d Vorticity Field: SP = 1.00. 
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Figure B.25a Constant Density Contours: SP = 1.00. 
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Figure B.25b .Velocity Field: SP = 1.00. 
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Figure B.25c Density Perturbation Contours: SP = 1.00. 
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Figure B.25d Voriicity Field: SP = 1.00. 
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Figure B.26a Constant Density Contours: SP = 1.00. 
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Figure B.26b Velocity Field; SP = l.OO. 
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Figure B.26c Density Perturbation Contours: SP = 1.00. 
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Figure B.26d Vorticity Field: SP = 1.00. 
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Figure B.27a Constant Density Contours: SP = 1.00. 
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Figure B.27b Velocity Field: SP = 1.00. 
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Figure B.27c Density Perturbation Contours: SP = 1.00. 
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Figure B.27d Vorticity Field: SP = 1.00. 
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Figure B.2Sa Constant Density Contours: SP = 1.00. 
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Figure B.28b Velocity Field; SP = 1.00. 
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Figure B.28c Density Perturbation Contours: SP = 1.00. 
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Figure B.28d Vorticity Field: SP = 1.00. 
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Figure B.29a Constant Density Contours: SP = 1.00. 
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Figure B.29b Velocity Field: SP = l.OO. 
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Figure B.29c Density Perturbation Contours: SP = 1.00. 
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Figure B.29d Vorlicity Field: SP = l.OO. 
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Figure B.30a Constant Density Contours: SP = 1.00. 
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Figure B.30b Velocity Field; SP = 1.00. 
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Figure B.30c Density Perturbation Contours: SP = 1.00. 
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Figure B.31a Constant Density Contours; SP = 1.00. 
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Figure B.31b Velocity Field; SP = 1.00. 
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Figure B. 31c Density Perturbation Contours: SP = 1. 00. 
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Figure B.3ld Vorticiiy Field: SP = l.OO. 
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T* VS H/BO (SP=1.00) 
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Figure B.32 . Experimental and Numerical Results for SP = 1.00. 
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Figure B.33a Constant Densitv Contours: SP = 0.50. 
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Figure B.33b Velocity Field; SP = 0.50. 
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Figure B.33c Density Perturbation Contours: SP = 0.50. 
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Figure B.33d Vorticicy Field; SP = 0.50. 
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Figure B.34a Constant Density Contours; SP = 0.50. 
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Figure B.34b Velocity Field; SP = 0.50. 
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Figure B.34c Density Perturbation Contours; SP = 0.50. 
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Figure B.34d Vorticity Field; SP = 0.50. 
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Figure B.35a Constant Density Contours: SP = 0.50.. 
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Tigure B.35b Velocity Field; SP = 0.50. 
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Figure B.35c Density Perturbation Contours: SP = 0.50. 



165 



Y/BO 



1 . 19040 





-5 -4 -3 -2 -1 0 



2 3 4 5 



X/BO 



Figure B.35d Vorticity Field; SP = 0.50. 
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Figure B.36a Constant Density Contours: SP = 0.50. 
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Figure B.36b Velocity Field: SP = 0.50. 
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Figure B.36c Density Perturbation Contours: SP = 0.50. 
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Figure B.36d 



Vorticity. Field; SP = 0.50. 
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Figure B.37a Constant Density Contours: SP = 0.50. 
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Figure B.37b Velocity Field: SP = 0.50. 
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Figure B.37c Density Perturbation Contours: SP = 0.50. 
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Figure B.37d Vorticity Field: SP = 0.50. 
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Figure B.3Sa Constant Density Contours: SP = 0.50. 
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Figure B.38b Velocity Field; SP = 0.50. 
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Figure B.3Sc Density Perturbation Contours: SP = 0.50. 
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Figure B.38d Vorticily Field: SP — 0.50. 



178 



3:57120 





-4 -3 -2 -1 0 1 2 3 4 5 

X/BO 

Figure B.39a Constant Density Contours: SP = 0.50. 
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Figure B.39b Velocity Field: SP = 0.50. 
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Figure B.39c Density Perturbation Contours; SP = 0.50. 
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Figure B.39d Vorticity Field: SP = 0.50. 
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Figure B.40a Constant Density Contours: SP = 0.50. 
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Figure B.40b Velocity Field: SP = 0.50. 
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Figure B.40c Density Perturbation Contours: SP 



0.50. 
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Figure B.40d Vorricicy Field: SP = 0.50. 
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Figure B.4la Constant Density Contours: SP = 0.50. 
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Figure B.41b Velocity Field; SP = 0.50. 
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Figure B.41c Density Perturbation Contours: SP = 0.50. 
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Figure B.41d Vorticity Field: SP = 0.50. 
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Figure B.42a Constant Density Contours: SP = 0.50. 
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Figure B.42b Velocity Field: SP = 0.50. 
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Figure B.42c Density Perturbation Contours: SP = 0.50. 
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Figure B.42d Vorticity Field; SP = 0.50. 
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Figure B.43a Constant Density Contours: SP = 0.50. 
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Figure B.43b Velocity Fieid; SP = 0.50. 
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Figure B.43c Density Perturbation Contours; SP = 0.50. 
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Figure B.43d Vorlicity Field; SP = 0.50. 
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Figure B.44 Experimental and Numerical Results for Various Initial Core Radii. 
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I'igure B.45a Constant Density Contours: SB = 0.50. 
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Figure B.45b Velocity Field; SP = 0.50. 
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Figure B.45c Density Perturbation Contours: SP = 0.50. 
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Figure B.45d Voriicity Field: SP = 0.50. 
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Figure 'B. 46a Constant Density Contours: SP = 0.50. 
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Figure B.46b Velocity Field: SP = 0.50. 
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Figure B.46c Density Perturbation Contours; SP 
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Figure B.46d Vorticiry Field: SP = 0.50. 
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Figure B.47a Constant Density Contours: SP = 0.50. 
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Figure B.47b Velocity F'ield; SP = O.fO. 
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Figure B.47c Density Perturbation Contours: SP = 0.50. 
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Figure B.47d Vorticity Field: SP = 0.50. 



211 



Y/BO 



1.78560 



T 



K 




-5 -4 -3 -2 -1 0 



2 3 4 5 



X/BO 



Figure B.4Sa Constant Density Contours; SP = 0.50. 
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Figure B.4Sb Velocity Field: SP = 0.50. 
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Figure B.48c Density Perturbation Contours: SP = 0.50. 
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Figure B.4Sd Voriicity Field: SP = 0.50. 
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Figure B.49a Constant Density Contours: SP = 0.50. 
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Figure B,49b Velocity Field; SP = 0.50. 
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Figure B.49c Density Perturbation Contours; SP = 0'.50. 
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Figure B.49d Vorticity Field: SP = 0.50. 
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Figure B.50a Constant Density Contours; SP = 0.50. 
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Figure B.50b Velocity Field: SP = 0.50. 



,221 



Y/BO 



T“ = 2.97600 




Figure BJOc Density Perturbation Contours: SP = 0.50. 
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Figure B.50d Vorticity Field: SP = 0.50. 
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Figure B.51a Constant Density Contours: SP = 0.50. 
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Figure B.51b Velocity Field: SP = 0.50. 
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Figure B.51c Density Perturbation Contours: SP = 0.50. 
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Figure B.51d Vorticity Field: SP = 0.50. 
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Figure B.52a Constant Density Contours: SP = 0.50. 
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Figure B.52b Velocity Field: SP = 0.50. 
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Figure B.52c Density Perturbation Contours: SP = 0.50. 
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Figure B.52d Vorticity Field; SP = 0.50. 
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